Three-dimensional gravity modelling of a Quaternary overdeepening fill in the Bern area of Switzerland discloses two stages of glacial carving

The geometry of glacial overdeepenings on the Swiss Plateau close to Bern was inferred through a combination of gravity data with a 3D gravity modelling software. The target overdeepenings have depths between 155 and > 270 m and widths between 860 and 2400 m. The models show incisions characterized by U-shaped cross-sectional geometries and steep to over-steepened lateral flanks. Existing stratigraphic data reveals that the overdeepenings were formed and then filled during at least two glacial stages, which occurred during the Last Glacial Maximum (LGM) within the Marine Isotope Stage (MIS) 2, and possibly MIS 6 or before. The U-shaped cross-sectional geometries point towards glacial erosion as the main driver for the shaping of the overdeepenings. The combination of the geometries with stratigraphic data suggests that the MIS 6 (or older) glaciers deeply carved the bedrock, whereas the LGM ice sheet only widened the existing valleys but did not further deepen them. We relate this pattern to the different ice thicknesses, where a thicker MIS 6 ice was likely more powerful for wearing down the bedrock than a thinner LGM glacier. Gravity data in combination with forward modelling thus offers robust information on the development of a landscape formed through glaciers.


S2: Estimation of the density of the Upper Marine Molasse bedrock using the Nettleton method
: Nettleton profile for the Belpberg. The trend line of the local gravity anomaly for the bedrock density of 2500 kg/m 3 parallels the regional gravity field trend. This suggests that the density value, for which there is the least correlation with the topography, is the best estimate for the bedrock density. The westernmost station is 1015 and the easternmost station is 3005 (Table S1) Because the Belpberg comprises a homogeneous suite of UMM sediments 1,2 (see main text), we designed our elevation profile ( Figure S1) such as that it is oriented as much perpendicular as possible to the structure of the Belpberg while remaining parallel to the stations. We assigned densities between 2300 kg/m 3 and 2670 kg/cm 3

to the Upper Marine
Molasse bedrock, the latter of which is the standard density commonly used in gravity research 3 . We calculated the resulting Bouguer anomalies for the target section across the Belpberg and iteratively increased the bedrock densities with 100 kg/m 3  -2 km to 20 km: 50 m DEM as mass lines.
We proceeded with these calculations until we found a best fit between the trend of the Bouguer anomalies across the Belpberg and the regional trend ( Figure S2). In particular, the Bouguer anomalies that are based on a bedrock density of 2670 kg/m 3 yields an inverse shape of the surface topography, suggesting that the density of the Upper Marine Molasse bedrock is lower than that. In contrast, the consideration of a density of 2400 kg/m 3 yields a positive correlation, implying that this density yields in an underestimation of the gravity contribution of the Upper Marine Molasse bedrock. We finally found a best fit between the trend of the Bouguer anomalies across the Belpberg and the regional trend if we consider a density of 2500 kg/m 3 for the Upper Marine Molasse bedrock ( Figure   S2). Note that there is a very local anomaly between points MP2001, MP2003 and MP2002 where we estimate a local density of 2450 kg/m 3 as indicated by the concave and convex triplet anomalies for density values of 2400 kg/m 3 and 2500 kg/m 3 , respectively. This suggests the occurrence of missing mass in a relative sense, which we attribute to the effect of moraine ridges in this region 1,2 .

S3: Development of the gravity modelling software package PRISMA
We developed a routine (see Python code in Supplementary Data referred to as Prisma) to address two major goals including (i) the assessment of the bulk densities of the Quaternary infill and (ii) the determination of the cross-sectional geometry of the overdeepenings based on 3D a-priori information that is already available ( Figure S4a To achieve these goals, we developed a code referred to as Prisma ( Figure S3 and Prisma in Supplementary Data). Prisma calculates the gravity effects of right-angled prisms with a given density on a station. It is thus used to calculate the gravity effect of simple structures, quarries and mountain ranges characterized by linear or rectangular geometries. We employ Prisma to calculate the topography correction for the Bouguer anomaly by representing the DEM cells with prisms. This routine can further be used to determine the cross-sectional geometry of the flanks of the bedrock beneath the surface, which is the scope of this paper. Prisma was also used to integrate the model and survey results into a GIS environment. Finally, Prisma serves to illustrate the first -order effects of a complex structure on a station. Please note that although we focus our work on reconstructing the cross-sectional and hence the 2D geometry of the target overdeepening, we calculated the 3D gravity effect of the overdeepening fill on our cross-section, because the gravity at a point is the consequence of the entire (and thus 3D) sedimentary architecture of the overdeepening fill in close proximity to our target cross-section.
We approximated the overdeepening fill with prisms, which are constrained using a-priori information that is based on the bedrock map underneath our target overdeepenings (see Reber and Schlunegger 4 for thickness maps). Please also note that other gravity toolboxes exist 6 that can be employed for various settings in Earth sciences. However, since Prisma is tailored towards addressing specific questions that are related to overdeepenings 3 (i.e., estimation of gravity contrasts between the bedrock and the overdeepening fill, estimation of the maximum depth of overdeepenings etc.), we employed our own code. Please see the main text for further explanations.  Figure S4a: Bedrock map 4 , which was used as a basis to allocate the prisms' positions in the Gürbe and Aare valleys. Please note that the top surface of the Aare valley prism is located at an elevation of 520 m a.s.l. , which corresponds to the station with the lowest elevation. Accordingly, for the Gürbe valley, we selected an elevation of 523 m a.s.l. The bedrock elevation model is openly accessible (SwissAlti3D, © swisstopo).

S4: Estimation of density contrast between the bedrock and the overdeepening fill and maximum thickness of overdeepening fill using Prisma
The estimation of the density of the Quaternary fill was accomplished indirectly 1 in the sense that the application of Prisma yields a density contrast between the bedrock and the sedimentary fill. For this purpose, for the Gürbe valley, we selected a prism with a width that fits both flanks and that runs parallel to the overdeepening as inferred by the bedrock topography map 4 ( Figure S4a). The top of the prism is located at an elevation of 523 m. a.s.l., the lowest elevation close to our gravity section, and the centre of the prism is slightly to the north of it to better comply with the slightly curved geometry of the Gürbe valley. Likewise, the 2 km-long distance of the Prisma to the North and to the South from the centre, respectively, is defined by the geometry of the western flank of the Gürbe valley. The width of 860 m from the centre is constrained by the flank of the valley. The elevation profile ( Figure S1) was positioned to be perpendicular to the prism and parallel as much as possible to the gravity stations. The 5° clockwise rotation of the Swiss coordinate system axes, along the centre of the prism, was done to comply with the conditions that the model requires right-angled prisms that are oriented parallel to both coordinates' axes of the prism 7,8 .
For the Aare valley, the prism was defined with its centre located at the same position as MP 3011, which corresponds to the valley centre. Both the NE and SW borders of the Aare valley prism were defined using the bedrock model from Reber and Schlunegger 4 as a criterion ( Figure S4a). In particular, we allocated the borders along the margins where the upper part of the slope of the overdeepening ends ( Figure S4a). The prism has a length of 2000 m in both directions and a total width of 2400 m. The elevation line is shifted to the North of the prism centre as it is the best compromise to fit the surface topography , which shows that the western flank of the Belpberg has a concave geometry. This profile also includes the stations at the centre of the valley. The prism is also offset to this line to best fit the information from the bedrock topography model of Reber and Schlunegger 4 .
A counterclockwise rotation of 33° was operated onto the Swiss coordinate system axes for the same reasons as the Gürbe valley.
The results show that for the Gürbe valley, a maximum residual anomaly of -2.9 mGal (see Figure 4 main text) is reached for a density contrast between 400 and 500 kg/m 3 , which correspond to densities of 2100 and 2000 g/cm 3 for the Quaternary sediments ( Figure S4b). All other scenarios are associated with thicknesses and densities for the Quaternary material that are beyond what has been estimated so far 5 Figure S4c). We note, however, that this thickness must be considered as a lower bound because available drilling data (see Figure 2c main text) imply a thicker Quaternary suite. Moreover, because it is a simplistic single prism model, the gravity effect on the flanks is overestimated, as they reach deeper than what drillings have disclosed, which further biases the estimation of the thickness of the Quaternary fill. Figure S4b: Plot of the Prisma model results that are used to determine the maximum Quaternary sediment thickness and the density contrast between the Quaternary sediments and the bedrock in the Gürbe valley, following Kissling and Schwendener 3 . Figure S4c: Plot of the Prisma model results that are used to determine the maximum Quaternary sediment thickness and the density contrast between the Quaternary sediments and the bedrock in the Aare valley, following Kissling and Schwendener 3 .

S5: Modelling of the residual gravity anomaly across the Gürbe valley using Prisma
The figures S5 present the model results of the residual gravity anomalies across the Gürbe valley. We started with a single prism ( Figure S5a), and the first model ( Figure   S5b) represents an approximation of the overdeepening fill, which can be described, in the simplest case, by one prism. The middle plot ( Figure S5c Figure 6 see main text). The DEM is openly accessible (SwissAlti3D, © swisstopo). Figure S5b: Modelling results where the overdeepening fill is approximated by one single prism. There is a large misfit between observed and modelled residual Bouguer anomalies, indicating that we considered too much Quaternary material on the overdeepening's flank, as the half width of the anomaly does not fit with a too high amplitude. The maximum anomaly, however, is properly modelled even for the stations where the anomaly is 90% of the maximum, indicating a likely wide geometry at depth for the overdeepening. The blue and green dots represent the observed gravity at the stations (see map above for legend) while the orange dots represent the modelled residual anomaly for each station. The blue line and blue dots show the elevation profile with the elevation of the stations, and the red dashed line and red diamonds indicate the bedrock model profile 1 and the bedrock depth from drillings, respectively.  give a satisfactory first approximation of the residual anomaly as the model anomaly has a higher amplitude than the overall observations ( Figure S6b). For this reason, we used two prisms. A wider top prism is used to approximate the shallower part of the overdeepening, whereas a narrower but much thicker prism is assigned for the depth of the overdeepening ( Figure S6c). We then tested the cases where we considered both V- dots. The drillings that reached the bedrock and that are used as constraints are marked in red, with the two yellow and white stars indicating two drillings that did not reach the bedrock but that are still used for constraining (Figure 3 main text) the geometry. Figure S6a: Map of the Aare valley prism. The prism is in orange, the blue dots show the stations. Because of the nonlinear organization of the stations, we do not distinguish between stations that are closer or farther away from the gravity profile. The red diamonds represent the drillings used to constrain the depth for the models. The two stars indicate drillings that provided a minimum depth constraint, and they are used for the sedimentary log of this paper ( Figure 6 see main text). The DEM is openly accessible (SwissAlti3D, © swisstopo).  Figure S6c: Results of the Aare valley modelled with two prisms (see Figure S6b for the legend). The results for this second model are much better, but the amplitudes for the maximum anomaly and the half-width are still too large. Yet it is much closer than the solution illustrated in F igure S6b (above). The flanks of the overdeepening have now been properly approximated, confirming the presence of a much narrower overdeepening geometry at depth than at the surface. (See Figure S6b for the legend). Figure S6d: Results of the V-shaped geometry model of the Aare overdeepening. The results show a fit of the general amplitude of the residual anomaly, however the fit around the maximum anomaly is not good, as the model results appear to underestimate the observed effect. The residual anomaly is slightly underestimated in this model. In the same sense, the anomaly on the western border is not properly modelled with this symmetric geometry. This model indicates that the slope of the Aare overdeepening is not as steep as the Gürbe valley overdeepening. It also indicates the need for an asymmetric geometry to approximate the anomaly pattern on the western flank. The input and output of this model is available in the supplementary dataset Prisma, and it is used as an examp le to run the code, labelled example2_vshape. (See Figure S6b for the legend). Figure S6e: Results of the U-shaped geometry model of the Aare overdeepening. Like above, the results here show a good fit of the general amplitude of the residual anomaly. In particular, the fit is reasonably good around the maximum of the anomaly, but it is worse on the western flank where the modelled amplitudes are too high. The effect on the eastern flank is slightly overestimated. This model suggests that a wide and flat geometry at the bottom of the overdeepening is very likely, and that the geometry of the eastern flank lies in between this one and the solution illustrated in Figure S6d (see above). It also shows the need for an asymmetric western flank, with an extended plateau at depth. (See Figure S6b for the legend).
Please note we are not capable of accounting for a possible depth dependency of the gravity values at least for the Holocene mud. However, the consideration of a lower density for the Holocene material and a higher density for the pre-LGM deposits will not alter our conclusion that the upper part of the target overdeepenings is much wider than the lower one (see Figures S6). Figure S4c suggest that a density contrast of 450 kg/m 3 could also be suitable to explain the maximum residual anomaly measured across the Aare valley. Using such a density contrast, however, does not yield a geometry that fits the both the gravity and the depth contrast offered by drillings on the eastern flank ( Figure S7). These results support our choice of the same density contrast for both valleys. Figure S7a: Similar to Figure S4c, this plot indicates the density contrasts and prism thickness needed to fit the maximum anomaly of the Aare overdeepening. Here we show that a contrast of 450 kg/ m3 could be a valid value to be used, because the prism thickness to explain the maximum residual anomaly would be similar to the proposed value in the bedrock model 4 . Figure S7b: Plot of the model that best fits the geometry of the Aare overdeepening on its eastern flank. This model is based on a density contrast of 450 kg/m3 as suggested by Figure S7a. Nevertheless, there is a misfit with the residual anomalies. In particular, this particular geometry does not provide an overdeepening volume that is large enough to match the observed gravity. Please see Figure S6 for the legend.